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c« ' Abstract 
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In this paper, we consider the problem of computing estimates of the domain-of- 
attraction for non-polynomial systems. A polynomial approximation technique, based 
on multivariate polynomial interpolation and error analysis for remaining functions, 

■ is applied to compute an uncertain polynomial system, whose set of trajectories con- 
tains that of the original non-polynomial system. Experiments on the benchmark non- 
polynomial systems show that our approach gives better estimates of the domain-of- 
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1. Introduction 



Stability for nonlinear control systems plays an important role in control system analysis and 
design. It will be very useful to know the domain of attraction (DO A) of an equilibrium point, 
however, this region is usually difficult to find and represent explicitly. Therefore, looking 
for underestimates of the DOA with simple shapes has been a fundamental issue in control 
system analysis since a long time. Among all the methods, those based on Lyapunov functions 
are dominant in literature [3, 4, 6, 7, 8, 11, 16, 19, 22, 20, 10, 21, 15]. These methods not 
only yield a Lyapunov function as a stability certificate, but also the corresponding sublevel 
sets as estimates of the DOA. 

For polynomial systems, many well-established techniques ([6, 8, 7, 11, 16, 19, 22, 20, 
10, 21, 15]) are available for computing estimates of DOAs. In [19], a method based on 
SOS decomposition was presented to find provable DOAs and attractive invariant sets for 
nonlinear polynomial systems. For odd polynomial systems, [6] employed an LMI-based 
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method to compute the optimal quadratic Lyapunov function for maximizing the volume 
of the largest estimate of the DOA. To obtain estimates of DOAs of uncertain polynomial 
systems, the authors of [8] used discretization (in time) to flow invariant sets backwards along 
the flow of the vector field. In [16], quantifier elimination (QE) method via QEPCAD was 
also applied to find Lyapunov functions for estimating the DOA. However, these methods 
cannot be applied directly in practice since most real systems are non-polynomial systems, 

1. e, their vector fields contain non-polynomial terms. For this kind of systems, only a few 
approaches have been proposed to deal with the DOA analysis. In [3, 4, 5], the author 
proposed an LMI technique through Taylor expansions as substitution for non-polynomial 
terms, and this technique can be generalized to compute estimates of DOAs for uncertain 
non-polynomial systems. In [23], an interval arithmetic approach was proposed. Recently, 
[18, 17] suggested a new method, based on quadratic Lyapunov function and the theorem of 
Ehlich and Zeller. 

In this paper, we will consider the problem of stability region analysis of uncertain non- 
polynomial systems. Through multivariate polynomial interpolation together with the inter- 
polation error analysis, we substitute a non-polynomial system as an uncertain polynomial 
system, whose set of trajectories contains that of the original non-polynomial system. By 
computing estimates of the DOA for the resulted uncertain polynomial system, we obtain 
estimates of the DOA for the original non-polynomial system. Our method is also applicable 
to the problem of searching for the largest possible underestimate of the DOA via a fixed 
Lyapunov function. Compared with the classical approximation by Taylor expansions, the 
error bound obtained using our suggested method is much sharper, which helps to yield a 
larger estimate of the DOA for a given non-polynomial system. 

The rest of the paper is organized as follows. In Section 2, some notions related to DOAs 
are presented. In Section 3, a polynomial approximation method, based on multivariate 
polynomial interpolation and interpolation error analysis, is proposed to substitute the non- 
polynomial functions as uncertain polynomials. In Section 4, bilinear SOS programming is 
applied to estimate DOAs of non-polynomial systems. In Section 5, experiments on some 
benchmarks are shown to illustrate our suggested method. Section 6 concludes the paper. 

2. Problem Formulation 

Consider an autonomous system 



where f : D C R™ — > W 1 is a continuous function defined on an open set D and f satisfies 
the Lipschitz condition: 



Denote by </>(£; x ) the solution of (1) with the given initial value x(0) = x . 

A vector x 6 I" is an equilibrium point of the system (1) if f(x) = 0. Since any 
equilibrium point can be shifted to the origin via a change of variables, we may assume 
without loss of generality that the equilibrium point of interest occurs at the origin. The 
equilibrium point of (1) is said to be stable, if for any e > there exists 6 such that 





f (x) — f (y) || < L||x — y || for all x, y 



e D. 
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Figure 1: Pendulum. 

whenever ||x || < S we have \\4>(t; x )|| < e for all t > 0; the point is said to be unstable if it 
is not stable; is asymptotically stable, if, in addition to being stable, there exists 6 such that 
lim i _ i . 00 (p(t; Xo) = whenever ||xo|| < 5; the equilibrium point is globally asymptotically 
stable, if, in addition to being stable, we have lim^oo (p(t; x ) = for all x G R n . 

Globally asymptotic stability is very desirable but is usually difficult to achieve. When 
the equilibrium point is asymptotically stable, we are interested in determining how far 
the trajectory of (1) can be from and still converge to as t approaches oo. This gives rise 
to the following definition. 

Definition 1 (Domain of Attraction). The domain of attraction (DO A) of the equilib- 
rium point for the system (1) is defined to be the set {x e W l \ lim^oo <p(t; x) = 0}. 

Usually, no algebraic description for DOAs is available. So researchers are mainly con- 
cerned with computing underestimates of the DOAs. Many well-established techniques 
([6, 8, 7, 11, 16, 19, 22, 20, 10, 21, 15]) are available for computing estimates of DOAs 
for polynomial (control) systems, i.e., autonomous systems with polynomial vector fields. 
However, in practice, many autonomous systems often contain non-polynomial terms in their 
vector fields. Below is an example. 

Example 1. [12, Example 1.2.1] Consider the simple pendulum shown in Figure 1. The 
motion of the pendulum is described by the following equation 

ml9 = —mg sin 9 — kW, 

where 9 denotes the angle subtended by the rod and the vertical axis through the pivot 
point, I the length of the rod, m the mass of the bob, g the acceleration due to gravity, and k 
the coefficient of friction. Let us take the state variables as x\ = 9, and x 2 = 9. Then the 
above equation is converted into a non-polynomial system 

r ±1 = x 2 , 

\ x 2 = -f sinxi - ^x 2 . 

□ 
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For the case of non-polynomial (control) systems, the problem of computing DOAs is still 
open, and only a few approaches have been proposed to deal with stability region analysis: in 
[3, 4, 5], the authors suggested a way to approximate non-polynomial vector fields by Taylor 
series expansion at the origin; in [23], an interval arithmetic approach for the estimation 
of the DOA was proposed; and recently, a method based on the theorem by Ehlich and 
Zeller was presented in [18, 17]. In this paper, we will apply polynomial approximation 
to transform a non-polynomial system into an uncertain polynomial system, whose set of 
trajectories contains that of the original non-polynomial system. Therefore, underestimate 
estimates of the DOA of the latter system yield those for the original non-polynomial system. 



3. Polynomial Approximation 

A key problem in estimating the DOA of a non-polynomial system is how to approximate the 
involved non-polynomial terms using polynomials, yielding an uncertain polynomial system 
with the equilibrium being kept. This problem is further reduced to the following problem. 



Problem 1. Let 0(x) : \1/ — > M be a non-polynomial function where f C ffi" is o bounded 
subset containing the origin 0. Given d G Z> 0; we will find a polynomial p(x) with de- 
gree d such that the error function r^(x) = 0(x) — p(x) satisfies r^(0) = and the value 
maXxg,!, is minimized. 

The classic method of polynomial approximation is Taylor expansions. Suppose <^(x) is a d 
times continuously differentiable in The Taylor expansion of 0(x) at the origin is 



\a\<d-l \p\=d H 



p(x) r d (x) 

for some £ G (0,x). In the above expression, p(x) is an approximate polynomial of 0(x) and 
the remainder term r^(x) is the error function of this approximation. Clearly, if the size of 
the region \1/ is small enough, the above Taylor expansion yields a tight bound of r<j(x) for all 
x£f. However, when the size of \l/ is large, the associated error bound may be too loose. 

To obtain a tighter bound, we will apply multivariate polynomial interpolation ([9]) to 
compute an approximate polynomial p(x) of <p(x) with a given degree d. Fix the graded 
lexicographic order in R[x]. For the function 0(x), one may find the minimal monomial x 7 
with 7 := ( 7l , . . . ,7n) G Z| , such that lim x ^ m J (0) ^ 0. Set V(x) = 0(x) ~ (/ ' (o) . Let d G 
Z> be such that d > \ j\ := ^27=1 7»- We construct a mesh Mon$ with mesh spacing s G M + 

and mesh points set x — {v x , v 2 , . . . , v fc } where k — ^ ^j- Like in [1], the meshes 

in our paper are either rectangular or simplicial. Then, we apply Lagrange interpolation to 
construct a polynomial p(x) as an approximation of ^(x) through the interpolation points Xi 
i.e., p(vi) = -0(vj), for i — 1, . . . , k. Next, we will compute a tight bound of the interpolation 
error function f(x) := ^(x) — p(x). Our idea is based on the following lemma. 
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Lemma 1. [26, Theorem 3] Let K C W 1 be a convex polyhedron, and vi, v 2 , . . . , and s 
be the vertices and diameter of K respectively. Suppose that ip : K — > R is a continuous and 
differential function on K, and A = sup xg ^ || V ¥>( x ) II ■ Then for all a\, a 2 , ■ ■ ■ , a& G swc/i 
£/jai ai + a2 + • • • + afe = 1, we have 

n 

|y?(x) - (ai^(vi) + a 2 <^(v 2 ) H h a fc ^(v fc ))| < — — -As. 

n + 1 

The following corollary gives an estimated bound of f(x) for x£f fl M. 

Corollary 1. Let s and x := { v i) v 2 5 • • • , v fc} ^ e ^ e J^es/i spacing and mesh points set 
of M, respectively. Suppose that p(x) is the interpolation polynomial of if)(x) through x, and 
f(x) = if}(x) — p(x) is the corresponding 6TTOT function. Let A — sviPxG^nA/ || \7 ^(-X-)II- Theu 

Tl 

|r(x) I < As for all x G M. 

n + 1 

Proof. Clearly, f(x) is a continuous and differential function on M, and 

f(vi) = f(v 2 ) = • • ■ = f(v fc ) = 0. 
Thus, according to Lemma 1, for all ai, a 2 , . . . , a v G R + such that a x + a 2 + • • • + a v = 1, 

|f(x) - (aif(vi) + a 2 r(v 2 ) H h a fc f(v fc ))| = |f(x)| < — — -As. 

n + 1 

□ 

Therefore, a non-polynomial function 0(x) can be relaxed to an uncertain polynomial, as 
shown in the following theorem. 

Theorem 1. For a non-polynomial function 0(x) with x G /e£ x 7 wit/i 7 G Z> 5e 
£/ie minimal monomial such that lim^o ^ 0. Lei M be a mesh on ^ wift mes/i 

spacing s G M+ and i/ie mesh point set x — { v i ; v 2, • • • ,^k}, in which k = 

Suppose that p(x) is the interpolation polynomial of ip(x) '■— ^^~^ ^ at x with degree < 
d — \ f(x) is the associated interpolation error function, and A = sup x6$nM || y f(x)||. 
Then for each x G \& D M we have </>(x) = p(x) + r^(x), where 

TL 

phi.) = 0(0) + p(x) x 7 and r^(x) = bx 1 with \u\ < As. (2) 

n + 1 

Clearly, the bound of the error r^(x) in (2) depends on the mesh spacing s, which can yield 
a tighter bound. Furthermore, the bound of r^(x) in (2) will converge to zero if d — > 00. 

Example 2. Consider the function (f>(x) = cosx with x G \& = [—1.2,1.2]. We want to 
compute a polynomial approximation for cosx. Based on Theorem 1, we can obtain an 
uncertain polynomial with degree 6, where 

p(x) = 1 - 0.5x 2 + 0.0416525x 4 - 0.00134386x 6 , 
r 6 (x) = ux 2 , -0.0000336 < u < 0.0000336. □ 



n + d- |7| 
n 
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4. Computation of Domain of Attraction 

In this section, we will consider an uncertain non-polynomial system of the form: 

x = f (x, 9) for all 9 G 6 C R*, (3) 

where 9 denotes a vector of uncertainty. Assume that the equilibrium point of interest occurs 
at the origin 0, i.e, f(0,6>) = for all 9 G 0. Denote by 0(t;x o ,#) the solution of (3) for 
the initial value x(0) = x and the uncertainty 9. The Domain of Attraction (DOA) of the 
system (3) is defined as 

{x G R n \ lim cf)(t- x, 9) = for all 9 G 6}. 

t — s-oo 

Lemma 1 in [19] can be modified a bit to compute underestimates of the DOA for (3) through 
Lyapunov functions, as described in the following theorem. 

Theorem 2. [22, Proposition 2.1] If there exists a continuously differentiable function V : 
M" R such that 

n v := {x G lR n : V(x) < 1} is bounded, 
V(0) = 0, 

v(x) > o, Vx g n v \{o}, 

V(x) = f£ • f (x, 9) < 0, Vx G fiy\{0}, wee 
then Qy an invariant subset of the DOA. 

When the equilibrium is asymptotically stable, the set Q c is clearly an underestimate 
of the DOA since every trajectory starting in Q c remains in Q c and approaches as t — > oo. 
And, if the equilibrium is globally asymptotically stable then the DOA will be the whole 
space E n . 

To enlarge the estimate fly given in Theorem 2, [19] defined a variable sized region 

Pp = {x G R n : g(x) < /3} 

with g(x) a fixed and positive definite polynomial in K[x], for instance, (?(x) = Y^=i x h an< i 
maximize /3 subject to the constraint Pp C f2y and the constraints in Theorem 2. Thus, the 
problem of computing Qy can be transformed into the following problem: 

max /? 

s.t. fly is bounded, 
V(0) = 0, 

V(x) > 0, Vx G fiy\{0}, 

7(x) = g • f (x, 9) < 0, Vx G Qy\{0}, W G 
^(x)</3t=^(x)<l. 



(4) 



Suppose that the non-polynomial system (3) has the following form 

k 

±i = fi(x,9) = /io(x,0) +5^/y(x,0)0 o -(x) J i = l,...,n, (5) 

3=1 
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where : R" x 1 4 K are polynomials in x for j = 0, . . . , k, and <pij : M n — > WL are 
non-polynomial functions for j = 1, . . . , k. Using the polynomial approximation technique 
in Section 3, we can replace each non-polynomial term <pij{x.) by an uncertain polynomial 
Pij(x.)+UijX. li: > with the bound \uij\ < bij. This gives rise to the following uncertain polynomial 
system: 

k 

±i = /i(x, 9) = /io(x, 9) + J2 fij(*i 9 )(Pij ( x ) + ^ x7lJ ) ( 6 ) 

3=1 

where |ity| < for z = 1, ... ,n. It is not hard to prove that the set of all trajectories of 
the system (5) is a subset of that of the system (6), and, consequently, the DOA of (6) is 
actually a subset of that of the system (5). Furthermore, the tighter the bound b^ is, the 
closer the DOA of (6) is to the DOA of the original system. 

Next, we consider how to find an optimal estimate of the DOA for the uncertain system (5) 
through computing the Qy in (4). Remark that the constraint V"(x) = ^ ■ f (x, 9) < in 
(5) may involve non-polynomial terms due to the existence of 0jj(x)'s. In this situation, 
replacing the above constraint by V^(x) = ^ • f (x, 9) < 0, the problem (4) can be relaxed as 
follow. 

Theorem 3. Let f(x, 9) = (/i(x, 9), . . . , / n (x, 9)) T with /j(x, 9) given in (6). IfV(x.) is a 
solution of the following polynomial optimization 

max B 

V,0 



s.t.V(0) = 0, 

V(x) > Vx G 1V\{0}, 



dV 

s(x) < 8 h V{k) < 1, 



§£ • f(x,0) < OVx g n v \{0},W9 e ey Uij g {±b iJ }, 



(7) 



then fly := {x 6 K" : ^( x ) < 1} is ar *> invariant subset of the DOA for (6), and therefore an 
invariant subset of DOA for (5). 

Proof. By construction of /i(x, 0)'s, we have 

dV - 

Vx e n v \{o},vo g e, 3u i3 e : f(x) = — • f(x,0). 

Clearly, if the constraints in (7) are fulfilled, the conditions in (4) also hold. Therefore, Q v 
is certainly a subset of Qy in (4). □ 

Assume that O is a semialgebraic set. For simplicity, we suppose O = {9 G O : > 0}. 
By rewriting the third, fourth and fifth constraints into equivalent empty set conditions, the 
condition (7) is transformed as 



V(0) = 0, 
{x G W 1 
{x G R n 
{x G W n 



V(x) < l,x^0,y(x) <O} = 0, 

V(x) < l,x ^ 0, > 0, §£ • f(x) > 0} = 0,Vuy G {±^}, 

^(x)< / s,y(x)>i,y(x)^i} = 0. 
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As stated in [19], Stengle's Positivstellensatz [2] be applied directly to solve (8). However, 
from the computational point of view, it is more efficient to replace all the inequations in 
(8) by inequalities of the form / > or / < 0. This can be done by introducing constants 
5 G R+ and polynomials of the form [19] l(x) = E" =1 ejxf l , where G R+ and m is assumed 
to be even. For example, by using S± and /i(x), the second condition in (8) can be relaxed as 

{x G R" : F(x) - 1 < 0, V(x) - Zi(x) + 5i < 0} = 0. 

Therefore, the problem (8) can be transformed into the following feasibility problem: 

max 8 
s.t. V(0) = 0, 



{xG 
{xG 
{xG 



i < o,y(x) 

av 



h(x) + 6i < 0} 



V(x) 

V(x) - 1 < 0, • f (x, 9) - Z 2 (x) + <5 2 < 0, ^(0) > 0} 
<7(x)-/3<0,y(x)-l-5 3 >0} 



V«»j G {=t6y} 



Suppose that £[x] is the set of sum of squares (SOS) polynomials in R[x]. Since the 
constraints in the above problem involve no equations and inequations, only a special case 
of Stengle's Positivstellensatz is needed, as shown in the following corollary. 

Corollary 2. Let F = {/j}j=i,„.,r be a set of polynomials in R[x]. The semi-algebraic set 

{xGM n :/ l (x)>0,z = l,...,r} 

is empty if and only if there exist polynomials sq, s\, . . . , s; G S[x] such that 



i 

5=1 
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w/iere G {/f 1 • • • : G Z> } . 

Applying Corollary 2, and removing all the crossing products of the involved inequalities, 
we obtain the following relaxed problem: 



max 8 

v,p 

s.t.F(O) = 0, 

<7 (x) + ai(x)(l - V(x)) + a 2 (x)(-y(x) + /x(x) - 8 X ) = 0. 
A (x) + Ai(x)(l - V(x)) + A 2 (x)(f ■ f>(x) + Z 2 (x) - 5 2 ) 
p (x) + Pl (x)(/3 - <?(x)) + p 2 (x)(^(x) - 1 - 6 3 ) = 0, 



0, 



(9) 



with er t (x), A t (x), p t (x) G S[x], t = 0, . . . , 2 and for any mj G {±6^}. 

Suppose that the Lyapunov function V(x) to be computed is a polynomial of degree d and 
has the form V(x) = J2 a c a^ i ° i where c a G R, x a = re" 1 • • ■ and a = (ai, . . . , a n ) G Z> 
with J^" =1 ctj < d. The decision variables of the problem (9) are j3 and the coefficients of all 
the unknown polynomials occurred in (9), such as V(x), cr t (x), A t (x) and p t (x). Clearly, some 
nonlinear terms which are products of the undetermined coefficients will occur in (9), which 
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yields a non-convex bilinear matrix inequalities (BMI) problem. To solve BMI problems, 
either a Matlab package PENBMI solver [13], which combines the (exterior) penalty and 
(interior) barrier method with the augmented Lagrangian method, can be applied directly, 
or an iterative method can be applied by fixing {3 and the polynomials alternatively, which 
leads to a sequential convex LMI problem. The reader can refer to [25] for more details. 

Remark that, the above proposed method is also applicable to computing the largest 
possible estimate of the DOA for a non-polynomial system (1) at through a fixed Lyapunov 
function V(x). Let fi c := {x G I" : V( x ) < c}. We will compute 

tt c * ={xG R n : V{x) < c*}, 

where c* = sup{c G R : V(x) < 0, Vx G f2 c \ {0}}. Due to the existence of non-polynomial 
terms, c* cannot be computed explicitly. Instead, we will compute the lower and upper 
bounds of c* as follows. Replacing the involved non-polynomial <pij by uncertain polynomials, 
computing the lower bound q of c* can be relaxed as the following problem: 

Cd = max c 

c ' x A (10) 
s.t. g • f(x) < 0, Vx G fi c \{0}, Vuij G {±b lJ }, 

where d is the degree of the interpolation polynomials. Clearly, Cd will converge to c* when 
d tends to oo. Next, we will search for a tight upper bound v<i of c*. To achieve this, let us 
look for Vd such that, for each G {±6jj}, the constant semi-algebraic system 

dV - 

V(x) -v d < OAx^ A — -f(x) > (11) 

has real solutions, which implies that Q Vd is not an estimate of the DOA. Based on bisection, 
Vd can be computed by Maple packages Regular Chains, DISCOVERER [24] and RAGLib [14]. 

5. Experiments 

Let us present some examples of DOA analysis of non-polynomial systems. 
Example 3. [4, Example 1] Consider a non-polynomial system 

f ±i — —X\ + x 2 + 0.5(e Xl — 1), 

\ X 2 = —X\ — X 2 + X\X 2 + X\ COSXi. 

To estimate the DOA of this system, we need to approximate the occurred non-polynomial 
terms e Xl and cosxi by uncertain polynomials. Based on the technique in Section 3, we 
obtain 

' e Xl = 1 + (1.0000004 + Ui)xi + • • • + 0.0014482244x?, 

cosxi = 1 - (0.5 + u 2 )x\ + 0.041669352^ - 0.0013878601x?, 
< -0.6 <xi< 0.6, 

-3 x 10~ 6 < ui < 3 x 10~ 6 , 
-1.2 x 10^ 6 <u 2 < 1.2 x 10~ 6 , 
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(a) boundaries of il C6 with c§ = 0.321064 and (b) boundaries of the estimates for dcgV^ = 2 
V = and dcg V = 4 



Figure 2: Results of Example 4. 

and the associated uncertain polynomial system. 

We first consider a fixed Lyapunov function V(x\,x 2 ) = x\ + x\. For the given degree 
6 of interpolation polynomials, after solving the corresponding SOS programming (10), we 
obtain the lower bound Cq = 0.321064 of c*, which is an improvement over the results in [4] 
with the lower bound 0.3210. Furthermore, by solving the problem (11) we obtain a tight 
upper bound v Q = 0.3216 of c*. 

Next, we estimate the DOA with variable Lyapunov functions. Suppose g{xi,X2) = 
x\ + x\. For deg V = 2, solving the SOS programming (9) with BMI constraints yields 

V{x u x 2 ) = 0.56678683^ + 0.23598133xix 2 + 0.92086339^, 
13 = 1.0453916, 

which is an improvement over the result from [4] where f3 = 1.0404. Similarly, for degV = 4 
solving the SOS programming (9) with BMI constraints yields (3 = 1.4001306 and 

V(x 1 ,x 2 ) = 0.068693712x? + • • ■ + 0.27723966^ + 0.2167918^, 

which is an improvement over the result from [4] where f3 = 1.2769. Therefore, fiy is an 
estimate of the DOA of the given system. Figure 2 shows the results obtained with Lyapunov 
functions of degrees 2 and 4. □ 

Example 4. [4, Example 2] Consider a non-polynomial system 

r ±i = x 2 , 

\ x 2 = —0.2x2 + 0.81 smxi cosxi — sinxi. 

Using the technique in Section 3, we obtain the approximations of the non-polynomial terms 
sin x\ and cos x\ as follows 

' sin 2xi = (2 + u x )x 1 - 1.3333091x? + 0.26625372x? - 0.023889715^, 

sinxi = (1 + u 2 )x 1 - 0.16666643x? + 0.008331760xf - 0.00019471928^, 
< -0.84 < xi < 0.84, 

-5.4 • 10~ 5 < ui < 5.4 • 10~ 5 , 
-1 • 10~ 7 < u 2 < 1 ■ 10~ 7 , 
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(a) boundaries of il C7 with C7 = 0.69922 and (b) boundaries of the estimates for dcg V = 2 
V = and dcg V = 4 



Figure 3: Results of Example 5. 

and the associated uncertain polynomial system. 

We first fix the Lyapunov function V{x\,x 2 ) — x\ + x\x 2 + ^x\. Let the degree of 
interpolation polynomials be 7. Solving the corresponding SOS programming (10), we obtain 
the results for the lower bound cj = 0.69922 of c*, which is an improvement over the result 
from [4] where the lower bound was 0.6990. Furthermore, by solving the problem (11) we 
obtain a tighter upper bound v 7 = 0.6998 of c*. 

We then estimate the DOA with variable Lyapunov functions. Suppose g(xi, x 2 ) = x\+x\. 
When degV = 2, by solving the SOS programming (9) with BMI constraints, we obtain 

V(x 1 ,x 2 ) = 1.01636667x? + 0.84993333xix 2 + 3.40233333:4 
P = 0.287706, 

which is an improvement over the result from [4] where j3 = 0.2809. Similarly, when deg^ = 
4, by solving the SOS programming (9) with BMI constraints, we obtain {3 = 1.92156 and 

V(xi, x 2 ) = 0.053849691x? + ■ ■ ■ + 0.11144243:4 + 0.058855229^, 

which is an improvement over the result from [4] where f3 = 1.1236. Therefore, Qy is an 
estimate of the DOA of the given system. Figure 3 shows the results obtained with Lyapunov 
functions of degrees 2 and 4. □ 

Example 5. [12, Example 2.2] Consider an uncertain non-polynomial system 

r xi = x 2 , 

\ x 2 = —6x2 — 10 sin £1, 

for 0.2 < 9 < 1. Based on the technique in Section 3, we obtain an approximation of the 
non-polynomial term sin X\ as follows 

sin 21 = (1 + ui)xi - 0.16664269012? + 0.008282118073xf - 0.0001751721223^, 
-2.4 < x 1 < 2.4, 

-4.365606 x 10~ 4 <u x < 4.365606 x 10~ 4 , 



and the associated uncertain polynomial system. 

Suppose g(x\, x^) = x\ + x\. For deg V = 4, solving the SOS programming (9) with BMI 
constraints yields (3 = 0.66552836 and 



V{x l ,x 2 ) = L1629845x? + ■ ■ • + 0.51014802x? + 0.010528:4 
Then Qy is an estimate of the DOA of the given system. □ 

6. Conclusion 

In this paper, we present a method on stability region analysis of non-polynomial systems via 
Lyapunov functions. A polynomial approximation technique, based on multivariate polyno- 
mial interpolation and error analysis, is applied to compute an uncertain polynomial system, 
whose set of trajectories contains that of the original non-polynomial system. To estimate 
DOA of the uncertain polynomial system, we apply Positivstellensatz to transform polyno- 
mial optimization problem into the corresponding (bilinear) sum of squares programming, 
which can be solved using the PENBMI solver or iterative method. Experiments on the 
benchmark non-polynomial systems show that our approach provides better estimates. 
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